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Preface 


The  purpose  of  this  study  was  to  evaluate  the  smooth 
particle  hydrodynamic  (SPH)  method  and  establish  its  credi¬ 
bility  as  a  tool  for  shock  wave  simulation.  SPH  is  a  relatively 
recent  numerical  method  and  shows  promise  for  2-  and 
3-dimensional  computations  of  gas-dynamical  flows. 

Evaluation  of  the  SPH  code  was  limited  to  one-dimensional 
shock  tube  problems  because  it  was  the  only  version  of  the 
code  available  for  study.  Therefore,  the  code  was  validated 
against  theoretical  predictions  of  shock  tube  behavior  and  a 
standard  Riemann  shock  tube  solution.  The  SPH  code  was  also 
compared  to  a  Lagrangian  hydrodynamic  code.  Por  similar 
resolution  of  shock  tube  density  features,  the  SPH  code  proved 
to  be  computationally  more  expensive.  It  also  displayed 
anomalous  behavior  at  the  contact  discontinuity.  This  was 
only  a  preliminary  evaluation  and  the  full  capability  of  the 
SPH  method  for  shock  wave  simulation  is  still  uncertain. 

I  gratefully  acknowledge  the  support  and  guidance  of  my 
faculty  advisor,  LCDR  Kirk  Mathews.  He  demonstrated  unexpected 
patience  and  understanding  during  the  low  points  of  my  research 
effort . 

Luke 
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Abstract 


A  smooth  particle  hydrodynamic  code  (SPHC)  is  evaluated 
for  performing  shock  wave  simulations  by  application  to  one¬ 
dimensional  shock  tube  problems.  Results  of  a  shock  tube  test 
case  with  a  compression  ratio  of  10  are  compared  against  a 
standard  Riemann  shock  tube  problem  and  theoretical  predictions 
of  shock  tube  behavior  to  validate  the  SPH  code.  A  Lagrangian 
hydrodynamic  code  is  validated  in  a  similar  fashion.  The 
resolution  capabilities  of  both  codes  are  compared  using  100, 
200  and  500  particles  for  SPHC  and  100,  400  and  800  cells  for 
the  Lagrangian  code.  The  SPH  code  exhibits  a  sharp  spike  in 
density  at  the  contact  discontinuity  for  a  shock  tube  test 
case  run  with  a  compression  ratio  of  100.  This  behavior  is 
not  reported  in  the  literature  and  not  seen  in  the  Lagrangian 
code  results.  Run  time  scaling  is  investigated  for  both  codes. 
SPHC  is  found  to  scale  between  NlogN  and  A/2,  where  N  is  the 
number  of  particles.  The  Lagrangian  code  scales  0(  N2). 
Computation  times  for  the  SPH  code  are  greater  than  run  times 
for  the  Lagrangian  code  by  a  factor  of  four  for  N  £  500  to 
achieve  similar  resolution  of  shock  tube  features. 
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SHOCK  TUBE  SIMULATION  BY  THE 
SMOOTH  PARTICLE  HYDRODYNAMIC 
(SPH)  METHOD 

I .  Introduction 

Approximately  50  percent  of  the  energy  released  by  a  fission 
weapon  in  air  at  moderate  altitudes  (40,000  ft)  or  lower 
contributes  to  production  of  a  shock  (or  blast)  wave  (1:7). 
Most  of  the  material  damage  resulting  from  the  nuclear  explosion 
is  attributable  to  the  blast  wave.  Blast  wave  phenomena, 
therefore,  are  of  great  interest  in  nuclear  weapon  effects  and 
nuclear  survivability  studies.  However,  gaining  experimental 
data  of  blast  wave  interaction  with  structures  is  extremely 
difficult.  All  aboveground  testing  of  nuclear  weapons  in  the 
United  States  ended  with  the  nuclear  weapons  test-ban  agreement 
of  1963  (2:273).  Alternative  methods,  therefore,  were 

developed.  Air  blast  testing  can  be  conducted  in  shock  tubes, 
simulated  using  large  high  explosive  (HE)  charges,  and  performed 
computationally.  The  high  cost  of  performing  shock  tube  and 
HE  testing  limits  the  frequency  and  variety  of  that  type  of 
testing . 

On  the  other  hand,  numerical  simulation  is  a  relatively 
inexpensive  method  for  studying  blast  wave  propagation  in  air 
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and  blast  loading  of  structures.  Numerical  hydrodynamic  methods 
offer  the  benefits  of  low  cost,  unlimited  parametric  studies, 
flexibility  of  boundary  configurations  and  conditions,  and 
repeatability  of  results.  Computational  methods  are  not  without 
their  problems,  though.  Computation  times  increase  dramati¬ 
cally  for  two-  and  three-  dimensional  shock  problems.  Tra¬ 
ditional  finite  difference  methods  involving  a  spatial  mesh 
suffer  from  mesh  tangling  or  inaccuracies  associated  with  the 
severe  distortion  of  the  mesh  for  multidimensional  problems. 
Eulerian  codes  may  require  such  fine  meshing  and  short  time 
steps  to  achieve  desired  resolution  in  the  solution  that 
computation  time,  and  therefore  cost,  becomes  prohibitive. 
Lastly,  numerical  solutions  are  only  as  good  as  the  physics 
that  goes  into  their  formulation.  Thus,  new  numerical  methods 
that  offer  performance  enhancements  and  increased  accuracy 
over  previous  hydrodynamic  codes  are  of  great  interest  and 
merit  evaluation.  One  of  these  new  methods  is  smooth  particle 
hydrodynamics  (SPH). 

1 . 1  SPH  Background 

Smooth  particle  hydrodynamics  was  developed  for  astro¬ 
physics  applications  involving  fluid  masses  moving  arbitrarily 
in  three  dimensions  in  the  absence  of  boundaries.  Lucy  (3) 
first  applied  SPH  to  the  problem  of  protostellar  fission. 
Gingold  and  Monaghan  (4)  extended  this  work  and  also  applied 
SPH  to  one-dimensional  shock  tube  problems  by  incorporating 
an  artificial  viscosity  into  the  equations  of  motion. 
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Experiments  showed  that  this  artificial  viscosity  term  produced 
negligible  oscillation  and  good  resolution  of  the  shock  front 
and  contact  discontinuity  (4:375).  Recently,  Benz  has  applied 
this  technique  to  planetary  (5)  and  stellar  (6)  collisions. 

Particle  methods  are  often  computationally  superior  to 
mesh-based  methods  when  computing  highly  asymmetric  three- 
dimensional  gas-dynamical  flows  (7:414).  Highly  distorted 
flows  and  large  voids  occur  during  impacts  or  mass  transfer 
in  stellar  collisions.  Lack  of  a  mesh  makes  SPH  particularly 
suited  for  simulating  such  problems  since  it  does  not  risk 
mesh  entanglement.  Also,  no  memory  or  computational  time  is 
wasted  by  having  a  large  number  of  empty  cells  just  in  case 
some  material  moves  into  them  (8:649)  as  in  Eulerian  schemes. 

Certain  limitations  of  these  methods  have  discouraged  their 
application  to  non-astrophysical  problems.  One  disadvantage 
of  SPH  is  run  time  scaling.  SPH  particles  interact  with  a 
variable  number  and  set  of  neighboring  particles  which  can 
lead  to  excessive  computation  (9).  Another  disadvantage  is 
particle  streaming  between  colliding  fluid  elements  (10).  This 
results  in  anomalous  fluid  interpenetration  which  does  not 
occur  in  a  real  physical  system.  Also,  for  fixed  mass  particles, 
low  density  regions  exhibit  poor  resolution  due  to  the  fewer 
particles  in  these  regions  relative  to  higher  density  ones. 
Lastly,  SPH  has  no  grid,  and  treatment  of  boundaries  has  been 
largely  ignored  until  recently  since  they  are  not  important 
in  astrophysical  applications. 
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Steilingwerf  (11)  reports  that  a  new  implementation  of  SPH 
by  the  Mission  Research  Corporation  (MRC)  overcomes  the  inherent 
disadvantages  of  SPH.  This  new  SPH  code,  SPHC,  has  been  tested 
by  MRC  against  analytic  solutions  and  other  hydrodynamic  codes 
on  rarefaction,  shock  tube,  blast  wave  and  collision  problems 
as  well  as  reproducing  laboratory  results  in  laser/target 
experiments.  It  is  this  code  which  will  be  evaluated  here. 

1 . 2  Problem  and  Scope 

The  objective  of  this  study  is  to  evaluate  the  merits  of 
the  SPH  method  as  a  simulation  tool  for  modeling  blast  wave 
phenomena  in  air  from  a  nuclear  explosion.  The  primary  goal 
is  to  evaluate  the  computational  performance  of  the  SPH  code 
in  resolution  capability  and  run-time  scaling  using  a  1-D  shock 
tube  test  case.  The  approach  is  to  compare  SPH  shock  tube 
results  to  a  conventional  (meshed)  1-D  Lagrangian  hydrodynamic 
code.  A  second  goal  is  to  apply  SPH  to  a  practical  application 
such  as  the  stagnation/amplification  of  a  shock  wave  in  a 
corner.  Due  to  difficulties  encountered  using  the  SPH  code, 
this  second  objective  was  not  met. 

1 • 3  Development 

A  theoretical  development  of  the  SPH  method  from  the 
integral  form  to  the  particle  approximation  is  presented  in 
Section  II.  Section  III  provides  details  of  the  shock  tube 
validation  for  the  SPH  and  1-D  Lagrangian  hydrodynamic  codes. 
A  comparison  of  the  resolution  capabilities  and  run  time  scaling 
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results  of  the  two  codes  is  detailed  in  Section  IV.  Section 
V  presents  an  analysis  of  the  SPH  code  such  as  its  unique 
features  and  problem  applications.  Conclusions  and  recom¬ 
mendations  are  presented  in  section  VI. 
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I I •  Theoretical  Development  of  SPH 

In  this  section  a  theoretical  development  of  the  SPH  method 
is  presented  based  on  the  references  by  Gingold  and  Monaghan 
(4,12). 

2 . 1  Theory 

SPH  is  a  free  Lagrangian  method  for  solving  the  conservation 
equations  of  hydrodynamics  by  replacing  the  continuum  of 
hydrodynamic  variables  with  a  finite  number  of  points.  SPH 
is  a  Lagrangian  method  in  the  sense  that  points  move  with  the 
fluid.  The  velocity  and  thermal  energy  at  these  points  are 
known  at  any  time.  A  mass  is  assigned  to  each  point,  so  they 
are  referred  to  as  particles.  Values  of  continuous  hydrodynamic 
variables  such  as  density  and  pressure  are  defined  by  an 
appropriate  average  over  the  distribution  of  particle  values 
using  a  special  weighting  function,  the  smoothing  function  (or 
interpolating  kernel). 

To  represent  a  continuous  flow  variable  at  r  by  a  smoothed 
local  average  over  a  distribution  of  particles,  consider  the 
kernel  estimate  of  the  function  /(r)  in  the  domain  D 

/«(r)  -  f  /(r')t/(r.r'./i)dr'  (1) 

J  D 

where 


6 


k  -  kernel  estimate 
r'  -  position  of  a  particle 
/(r')  -  the  function  interpolated 

W  -  the  interpolating  kernel 

The  interpolating  kernel  has  the  following  properties: 

(a)  jV(r.r',ft)dr'  ->  l  as  ft -» 0.  (2) 

D 

(b)  If  /(r)  is  a  continuous  function 

/«(D  -*  /(r)  as  ft  -*  0.  (3) 


Therefore,  the  interpolating  kernel  acts  like  a  delta  function, 
and  it  does  so  more  closely  as  ft-*  0.  The  smoothing  length  ft 
is  chosen  so  that  it  determines  the  extent  to  which  W  confines 
the  major  contribution  to  /K(r)  to  the  neighborhood  of  r-r' 
(12:423) . 


Although  there  are  many  possible  kernels  to  choose  from, 
only  two  kernels  are  commonly  used  (6:616):  the  Gaussian  and 
exponential  kernels.  The  Gaussian  kernel  is  given  by 

exp[-(r-r')2//i2] 


W(T,r',h) 


(n/tz) 


2^/2 


(4) 


whereas  the  exponential  kernel  is  given  by 


exp[-(|r-  r'|)/h] 
8  lift3 


(5) 


Assuming  that  we  have  a  set  of  N  points  r,.r2 . rN  dis¬ 

tributed  in  space  according  to  the  number  density  n(r),  /*(r) 
from  Eq  (1)  can  be  approximated  by 

/.( r)  -  (6) 

This  expression  for  /K(r)  is  the  general  SPH  interpolation 
formula.  The  usual  case  for  particle  methods  is  to  define  the 
mass  density  p(r)  as 

p(r)  -  m(r)n(r)  (7) 

where  m(r)  is  the  mass  of  each  particle.  For  most  purposes 
it  is  assumed  that  l/(r,ry,  /i)  ■  l/(|r  -  ry[  ,/i).  Eq  ( 6)  then  becomes 

/«(*■)  -  Z/z^I/dr-rJ./i)  (8) 

7-1  P7 

where,  for  any  function  B,  B,  ■  5(r;). 

The  gradient  of  a  function  is  given  by 

V/.(r)  -  £/y^ivi/(|r-r;|.h)  (9) 

7-1  P7 

The  next  step  in  developing  the  SPH  method  is  to  construct 
equations  suitable  for  numerical  work.  Each  of  the  exact 
equations  is  multiplied  by  the  interpolating  kernel  and  then 
integrated  over  the  domain  where  a  solution  is  required. 

Integration  by  parts  gives  the  equations  for  numerical  work. 

Monaghan  and  Gingold  (4)  recommend  simply  dropping  nonlinear 
boundary  terms  since  the  kernel  should  mimic  a  delta  function 
which  goes  to  zero  sufficiently  far  from  the  particles  rep- 
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resenting  the  edge  of  the  system.  However,  if  one  is  to 
consider  problems  where  the  particles  interact  with  a  boundary 
condition,  such  as  an  externally-applied  pressure,  one  must 
specifically  include  the  boundary  terms  (13:5).  The  SPH  code 
evaluated  in  this  study  includes  the  nonlinear  boundary  terms 
and  will  be  discussed  in  the  following  sections. 
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III.  1-D  Shock  Tube  Validation 


The  equations  governing  fluid  behavior,  theoretical  pre¬ 
dictions  of  shock  tube  behavior  and  validation  of  the  SPH  and 
Lagrangian  hydrodynamic  codes  are  presented  in  this  section. 
3 . 1  Governing  Equations 

The  basic  mathematical  model  of  the  physical  laws  governing 
the  motion  of  an  ideal  or  inviscid  fluid  is  the  set  of  Euler 
equations.  These  equations  of  motion  are  derived  by  applying 
the  principles  of  conservation  of  mass,  momentum  and  energy 
to  a  differential  volume  of  compressible  fluid. 

For  most  practical  purposes,  the  conservation  equations 
are  expressed  as  a  set  of  nonlinear  partial  differential 
equations  (PDEs).  The  resulting  three  equations  are  functions 
of  four  variables  which  describe  the  fluid  medium:  pressure 
p,  density  p,  fluid  velocity  u,  and  the  internal  energy  per 
unit  mass  e.  Adding  an  equation  of  state  of  the  f  orm  e  - /(p,p) 
completes  the  system  of  equations  and  allows  for  a  solution. 
This  system  of  equations  describes  the  variations  in  fluid 
properties  throughout  time  and  space.  Also,  (assuming  that 
pressures,  velocities  and  temperatures  change  sufficiently 
slowly  with  distance)  heat  conduction,  viscous  and  all  other 
internal  irreversible  processes  are  neglected  (14:6). 


10 


The  Euler  equations  in  conservation  form  for  one  dimension 


are : 


d_U  +  df(U) 
dt  dx 


(ID 


where 


~  p 

pu 

pu 

f(U)  - 

pu2+  p 

_p®<_ 

_pu(e,+  p/p) 

and  e,  is  the  total  energy  per  unit  mass  given  as 

u2 

e,  m  e  +  -  (12) 

Internal  energy  is  related  to  pressure  by  the  equation  of  state 
for  a  perfect  gas 

p  -  (v-l)pe  (13) 

Calculating  the  flow  properties  across  the  shock  front  requires 
applying  the  Rankine-Hugoniot  conditions  for  changes  across 
a  normal  shock. 

A  shock  tube  can  be  modelled  using  these  hyperbolic  con¬ 
servation  equations,  Eq  (11) ,  as  a  set  of  initial -value  problems 
with  piecewise  constant  initial  data 

(/(x  ,0)  m  U  L  ( for  x  <  0 ) ;  i/(x,0)-£/,  (/or  x>0)  (14) 

separated  by  a  discontinuity  at  x  -  0.  This  special  ini¬ 
tial-value  problem  is  defined  as  a  Riemann  problem  (15:357). 
The  shock  tube  problem  is  a  special  Riemann  problem  with 
uL  m  u#  ■  0  (16:664).  The  Riemann  shock  tube  problem  has  an 
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exact  solution  and  thus  provides  a  benchmark  for  approximate 
solutions.  Figure  1  compares  exact  and  approximate  solutions 
to  a  Riemann  problem  for  air  with  a  density  ratio  of  8:1. 


Figure  1.  Exact  and  Approximate  Solutions  to  a  Riemann 
Problem  (15:368) 


An  algorithm  for  calculating  approximate  solutions  to  a  Rie¬ 
mann  problem  was  presented  by  Roe  in  Reference  15.  However, 
the  equations  for  the  exact  solution  were  not  included;  the 
algorithm  for  the  exact  solution  was  taken  from  Sod  (17). 
Thus,  Figure  1  only  provides  a  qualitative  example  of  the 
expected  density  profile  in  a  shock  tube  for  time  t  >  0. 
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3.2 


Figure  2a  shows  a  simple  shock  tube.  A  gas-tight  diaphragm 
separates  the  shock  tube  into  two  sections.  One-half  of  the 
tube,  known  as  the  compression  chamber,  contains  the  compressed 
gas  at  a  pressure  p3,  which  is  in  excess  of  the  pressure  p0 
of  the  gas  in  the  other  half  of  the  tube  known  as  the  expansion 
chamber . 


Figure  2.  Shock  Tube,  a)  assembly,  b)  pressure  history  at 
t  =  0,  c)  t  >  0  (14:29) 


Upon  shattering  the  diaphragm,  a  shock  wave  propagates 
into  the  expansion  chamber  while  a  rarefaction  wave  propagates 
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back  into  the  compression  chamber.  Figures  2b  and  2c  show  the 
pressure  profile  along  the  shock  tube  before  and  after  the 
diaphragm  has  shattered.  The  dotted  line  in  Figure  2c  indicates 
the  position  occupied  by  that  gas  which  was  originally  at  the 
diaphragm.  The  gas  to  the  right  has  been  compressed  and  heated 
by  the  shock  wave  while  the  gas  to  the  left  of  this  line  has 
expanded  from  the  compression  chamber  and  has  therefore  been 
cooled.  This  position  is  a  contact  discontinuity.  At  this 
point  there  will  be  a  discontinuity  in  gas  temperature  and 
density.  Pressure,  however,  remains  constant  across  this 
discontinuity  as  seen  in  Figure  2c.  If  this  were  not  the  case, 
a  shock  front  would  form  at  this  location  as  well.  Velocity 
also  remains  constant  across  the  contact  discontinuity. 

When  the  shock  wave  encounters  the  rigid  wall  at  the  end 
of  the  shock  tube,  it  undergoes  a  reflection  and  will  be 
reflected  as  a  shock  wave.  Likewise  the  rarefaction  wave 
travelling  into  the  compression  chamber  also  encounters  a  rigid 
wall  and  will  be  reflected  as  a  rarefaction  wave. 

3 • 3  SPHC  Shock  Tube  Validation 

A  1-D  shock  tube  test  case  was  run  with  SPHC  to  provide 
confidence  in  its  output.  The  program  default  parameters  were 
used.  Table  1  lists  the  initial  conditions.  It  is  assumed 
that  the  gas  on  both  sides  of  the  diaphragm  are  in  thermal 
equilibrium  and  therefore  are  both  at  the  same  temperature. 
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The  compression  ratio  r]  is  defined  as  the  ratio  of  the  gas 
density  in  the  compression  chamber  to  the  gas  density  in  the 
expansion  chamber  at  time  =  0. 


Table  1 

SPHC  Test  Case  Conditions 


No.  Of  Particles 

100 

Gas 

Air 

Initial  Gas  Density 

1.004  kg/m3 

Compression  Ratio 

10 

Temperature 

300  K 

Shock  Tube  Length 

1.0  cm 

Diaphragm  Location 

0.4  cm 

Smoothing  Length 

0.015 

Art.  Vise.  Parameters 

a 

1.0 

P 

1.0 

The  results  of  the  SPHC  shock  tube  test  case  are  displayed 
in  Figures  3a  and  3b  for  t  =  0  and  t  =  lxlO'6  time  units 
respectively.  The  units  for  time  are  not  reported  by  the  code. 
Also,  SPHC  does  not  explicitly  calculate  pressure  values  for 
output  so  only  density  profiles  are  plotted.  In  Figure  3a  a 
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slight  deviation  from  a  vertical  slope  in  density  occurs  at 


the  diaphragm  location.  This  results  from  the  resolution 


capability  of  the  code  for  finite  numbers  of  particles. 


Increasing  the  number  of  particles  used  in  the  computation 


moves  the  density  profile  closer  to  the  expected  vertical  slope 


at  the  diaphragm  location. 


Figure  3b  displays  the  density  profile  at  lxlO’6  time  units 


atter  the  diaphragm  is  shattered.  It  exhibits  the  same 


qualitative  density  profile  as  that  seen  in  Figure  1.  However, 


the  changes  in  density  at  both  the  contact  discontinuity  and 


the  shock  front  are  smeared  and  extend  over  a  smoothing  length 


range  of  ~5 h  and  ~3 h  respectively.  Increasing  the  number  of 
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Figure  3.  a)  SPH  shock  tube  test,  100  particles, 
h  -  0.015,  t  =  0 
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Figure  3.  b)  SPH  shock  tube  test,  100  particles, 
h  -  0.015,  t  =  1  x  10'6 

particles  from  100  to  500  much  appears  to  improve  the  resolution 
as  seen  in  Figure  4.  Now,  the  changes  in  density  at  the  contact 
discontinuity  and  shock  front  are  quite  sharp  and  only  a  slight 
rounding  is  apparent  at  the  contact  discontinuity.  However, 
increasing  the  number  of  particles  by  a  factor  of  5  resulted 
in  a  reduction  of  /i  by  1/5.  The  density  is  still  smeared  over 
the  same  h  at  both  the  shock  front  and  contact  discontinuity. 
The  only  difference  is  that  the  smoothing  length  is  reduced. 
The  artificial  viscosity  parameters  must  be  varied  to  reduce 
smearing  of  the  shock  front  and  contact  discontinuity. 
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It  is  apparent  that  the  1-D  shock  tube  version  of  SPHC 
produces  results  that  qualitatively  agree  with  theoretical 
predictions  and  an  exact  solution  to  a  standard  Riemann  shock 
tube  problem.  Now  that  confidence  in  the  SPH  code  is  estab¬ 
lished,  shock  tube  results  of  the  Lagrangian  hydrodynamic  code 
will  be  analyzed  in  a  similar  fashion. 


Figure  4.  SPH  shock  tube  test,  500  Particles, 
h  -  0.003,  t  =  1  x  10'6 


3 . 4  Lagrangian  Shock  Tube  Val idation 

A  1-D  shock  tube  test  case  was  run  using  a  modified  version 
of  a  Lagrangian  hydrodynamic  algorithm  developed  by  Bridgman 
(18).  The  test  case  was  run  at  the  same  initial  conditions 
as  those  listed  in  Table  1  except  200  computational  cells  were 
used  instead  of  100  particles  and  the  artificial  viscosity 
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parameter  a  m  2*  Thermal  equilibrium  was  assumed  to  exist 
between  the  gases  on  both  sides  of  the  diaphragm.  This 
assumption  will  continue  to  be  applied  to  all  test  cases  in 
this  study.  One  advantage  of  the  Lagrangian  hydrocode  over 
SPHC  is  that  it  explicitly  calculates  values  of  pressure  for 
output . 

Figures  5  and  6  show  the  pressure  profiles  at  times  0  and 
6x  10’4 seconds ,  respectively,  after  the  diaphragm  is  shattered. 
These  results  almost  perfectly  model  the  predicted  pressure 
behavior  displayed  in  Figures  2b  and  2c.  There  is  a  slight 
dip  in  pressure  followed  by  some  minor  oscillations  at  the 
tail  of  the  rarefaction  wave  in  Figure  6.  This  behavior  has 
been  seen  in  other  numerical  methods  incorporating  an  artificial 
viscosity  term  (4:382).  The  density  profile  also  displays 
some  oscillatory  behavior  as  seen  in  Figure  7.  The  density 
changes  at  the  contact  discontinuity  and  shock  front  are  sharp 
and  almost  vertical  in  slope.  A  small  smearing  of  the  contact 
discontinuity  and  shock  front  will  occur  due  to  the  artificial 
viscosity . 

The  Lagrangian  hydrocode  showed  excellent  agreement  with 
predicted  behavior  for  both  density  and  pressure  profiles. 
Confidence  in  the  ability  of  both  codes  to  produce  credible 
results  for  shock  tube  problems  has  been  established. 


Figure  7,  Hydrocode  density  profile,  t  =  6xl0~4sec 
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IV .  Results 


A  compression  ratio,  T|,  of  10  is  not  typical  of  the  mag¬ 
nitudes  at  which  actual  shock  tube  tests  are  run.  The  Defense 
Nuclear  Agency  (DNA)  is  developing  a  Large  Blast/Thermal 
Simulator  (LB/TS)  with  a  driver  gas  at  a  pressure  of  1740  psig 
(19).  Ambient  pressure  is  14.7  psig,  so  this  equates  to 
r|  m  120.  The  following  shock  tube  test  cases  were  run  at 
T)  -  100  to  better  simulate  actual  test  conditions. 

Comparisons  between  SPHC  and  the  Lagrangian  hydrocode  were 
made  when  prominent  features  such  as  the  shock  front  reached 
the  same  spatial  position  for  both  codes.  Time  comparisons 
could  not  be  made  since  the  Lagrangian  hydrocode  had  times 
approximately  500  times  greater  than  SPHC.  For  example,  the 
shock  front  reached  x  =  0.8  cm  at  t  =  0.5  msec  for  the  hydrocode 
whereas  the  SPHC  code  reached  the  same  position  at  t  =  lxlO’6 
time  units.  SPHC  does  not  provide  the  units  for  time  in  the 
output  nor  did  the  users  manual  specify  the  units. 

4 . 1  Resolution  of  Density  Profiles 

Figures  8a,  b  and  c  show  the  density  profiles  for  SPHC 
using  100,  200  and  500  particles,  respectively.  As  expected 
the  resolution  of  the  density  change  at  both  the  contact 
discontinuity  and  shock  front  improve  with  increasing  number 
of  particles.  However,  an  anomalous  spike  occurs  at  the  contact 
discontinuity  that  was  not  seen  in  the  SPHC  results  at  rj  -  10. 
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Figure  8.  Density  profiles  for  SPH  at  r|  =  100  using 
a)  100  and  b)  200  particles 
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Figure  8.  c)  Density  profile  for  SPH  using  500  particles 

This  feature  was  not  reported  in  the  references  on  shock 
simulation  by  SPH  (4,10)  since  they  conducted  the  shock  tube 
tests  at  r|  -  4. 

The  spike  in  density  to  the  left  of  the  contact  discontinuity 
may  result  from  a  piling  up  of  particles  in  this  high  density 
region.  Stellingwerf  (11)  reports  that  SPH  codes  were  normally 
limited  to  a  density  contrast  of  about  a  factor  of  3  due  to 
fixed-size  particles.  The  density  contrast  across  the  contact 
discontinuity  is  a  factor  of  5.  SPHC  has  a  particle  divi¬ 
sion/recombination  option  which  is  supposed  to  maintain  res¬ 
olution  in  regions  of  high  density  contrast.  This  option  was 
not  used  in  this  test  case  due  to  time  limitations. 
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Figures  9a,  b  and  c  show  the  density  profiles  for  the 
Lagrangian  hydrocode  using  100,  400  and  800  cells,  respectively. 
The  hydrocode  does  a  poor  job  of  resolving  the  contact  dis¬ 
continuity  and  the  shock  front  using  100  cells  as  seen  in 
Figure  9a.  Increasing  the  number  of  cells  to  400  dramatically 
improves  the  resolution  although  a  small  dip  in  density  occurs 
at  the  tail  of  the  rarefaction  wave  as  seen  in  earlier  results. 
Another  increase  to  800  cells  results  in  a  slight  improvement 
in  the  density  dip.  As  in  the  previous  section,  an  artificial 
viscosity  parameter  of  a  ■  2  was  used. 
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4 . 2  Run  Time  Scaling 

Computation  times  for  individual  runs  for  both  SPHC  and 
the  Lagrangian  hydrocode  are  listed  in  Table  2. 


Table  2 

Run  Times  for  SPHC  &  Lagrangian  Hydrocode 


SPHC  a 

Lagrangian 

Hydrocode  6 

#  particles 

run  time 

#  cells 

run  time 

(min) 

(min) 

50 

3 

100 

1.3 

100 

9 

200 

5.1 

200 

22 

400 

20 

500 

93 

800 

82 

a  Sun  SPARC 2  workstation 
b  386SX,  16  MHz  PC 


The  500  particle  SPHC  run  and  the  400  cell  hydrocode  run 
give  qualitatively  comparable  resolution  of  the  density  pro¬ 
file.  However,  SPHC  requires  over  four  times  the  computation 
time.  The  references  on  SPH  methods  support  this  finding  and 
only  claim  that  SPH  schemes  are  computationally  superior  in 
2-  and  3-D  geometries  (7).  Computation  times  between  SPHC  and 
the  hydrocode  may  be  even  greater  since  SPHC  was  run  on  a  Sun 
SPARC 2  workstation  while  the  hydrocode  was  run  on  a  386-SX, 
16  MHz  PC.  The  Sun  is  an  estimated  5-10  times  faster  than 
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the  PC.  Therefore,  a  conservative  estimate  is  that  the 
Lagrangian  hydrocode  is  running  twenty  times  faster  than  SPHC 
for  comparable  resolution. 

Computation  time  for  the  hydrocode  increases  by  a  factor 
of  four  for  every  doubling  of  the  number  of  cells.  This  equates 
to  run  time  scaling  of  0(/V2).  The  literature  (9,  20)  predicts 
that  SPH  methods  using  a  hierarchical  tree-structured  technique 
for  force  calculations  between  N  particles  will  scale  as 
O(iVlogiV).  Figure  10  shows  that  the  SPHC  run  times  are 
N  log  N  <  t  <  N2. 

The  Lagrangian  hydrocode  proved  to  require  less  computation 
time  than  SPHC  for  N  £  500.  However,  N2  scaling  has  a  steeper 
slope  than  N\ogN  scaling,  so  for  problems  involving 
N  -  103- 104,  the  Lagrangian  hydrocode  becomes  computation¬ 

ally  more  expensive. 
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SPHC  run  time  scaling 
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V.  Analysis  of  SPHC 

A  discussion  of  the  unique  features  of  the  SPH  code  is 
presented  below.  Also  presented  are  the  problems  encountered 
in  evaluating  the  code  and  the  change  in  scope  of  the  study. 
5 . 1  Smooth  Particle  Hydrodynamic  Code  ( SPHC) 

Mission  Research  Corporation  (MRC)  of  Albuquerque,  NM  has 
developed  a  benchmark  SPH  code,  SPHC,  as  a  testing  arena  for 
the  SPH  technique.  What  sets  this  code  apart  from  earlier 
pioneering  efforts  in  SPH  is  its  emphasis  on  simulating  non- 
astrophysical  problems.  SPHC  incorporates  many  features  to 
overcome  the  deficiencies  found  in  earlier  codes  that  limited 
their  use  primarily  to  the  astrophysics  realm. 

5.1.1  SPHC  Features .  SPHC  addresses  the  problems  of  run 
time  scaling,  limited  density  range,  fluid  interpenetration 
and  treatment  of  boundaries  (11:2-4). 

1)  Run  time  scaling.  To  reduce  computational 
requirements,  SPHC  uses  a  hierarchical  tree  scheme 
for  neighbor  location  which  ensures  scaling  of 
A/logA/  instead  of  Nz  where  N  is  the  number  of  par¬ 
ticles. 

2)  Density  range.  Original  SPH  methods  used  fixed 
size  particles  which  limited  density  contrast  to  a 
factor  of  3.  MRC  uses  variable  size  particles  and  a 
new  technique  that  allows  particle  division  in  low 
density  regions  to  maintain  resolution. 
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3)  Interpenetration.  A  specialized  tensor  artifi¬ 
cial  viscosity  (10:591-595)  is  used  to  address  the 
problem  of  particle  streaming  at  interfaces.  This 
problem  is  especially  severe  in  shocks  and  colli¬ 
sions  . 

4)  Boundaries.  SPH  is  a  gridless  technique  and 
boundaries  cannot  be  identified  with  particular  par¬ 
ticles.  Therefore,  SPHC  incorporates  a  wall  bound¬ 
ary  based  on  terms  normally  dropped  in  the 
integration-by-parts  step  of  converting  equations  to 
numerical  format.  A  ghost  particle  boundary  is  also 
used  for  reflecting  and  periodic  cases. 

In  addition  to  the  above,  SPHC  has  several  other  new 
features.  It  is  written  in  C  and  can  be  run  on  Cray,  Sun  and 
even  MS-DOS  environments.  SPHC  runs  in  1-,  2-,  or  3-dimension 
modes  as  well  as  in  spherical  and  cylindrical  geometries.  MRC 
has  also  included  new  models  to  treat  electron  thermal  con¬ 
duction,  single  group  radiation  diffusion,  laser  deposition 
and  inertial  confinement  fusion  (ICF)  physics. 

5.2  Code  Evaluation. 

SPHC  is  reported  by  Stellingwerf  (11:4)  to  have  been  "thoroughly 
and  carefully  tested  against  analytic  solutions  and  other  codes 
on  rarefaction,  shock  tube,  blast  wave,  exponential  atmosphere, 
and  collisions  problems."  Two  factors  made  SPHC  particularly 
attractive  for  this  study:  1)  standard  test  problems  have 
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already  been  set  up  and  executed  for  2-  and  3-D  blast  wave 
calculations  and  2)  SPH  methods  are  often  computationally 
superior  to  mesh-based  methods  for  2-  and  3-D  applications. 

The  initial  objectives,  therefore,  of  this  study  were 
twofold.  First,  scope  the  computational  limits  of  this  particle 
method  and  evaluate  its  potential  as  a  simulation  tool  for 
studying  blast  wave  phenomena.  Second,  apply  the  code  to  the 
problem  of  shock  wave  stagnation/reflection  in  a  corner. 
Unfortunately,  the  initial  objectives  could  not  be  achieved 
due  to  unforeseen  problems  with  the  code.  This  will  be  discussed 
in  the  following  paragraphs. 

The  SPHC  production  code  is  not  a  general  purpose  code. 
Different  hydrodynamic  problems  such  as  rarefaction,  shock 
tube  or  blast  wave  cannot  be  executed  using  the  same  main 
module  code,  sphc.c.  Although  all  use  the  same  underlying 
physics,  MRC  has  developed  individual  SPHC  versions  for  each 
test  case.  However,  setup  menus  are  included  during  program 
initialization  to  allow  adjustment  of  parameters  such  as  number 
of  time  steps,  boundary  types,  physics  models  and  equation  of 
state  parameters.  Making  changes  to  the  test  problem  con¬ 
figuration  beyond  available  parameter  adjustments  requires 
modifying  the  module  setup  file,  sph_init.c,  and  recompiling 
the  code. 

The  SPHC  version  employed  in  this  study  is  a  version  of 
the  standard  Riemann  shock  tube  test  case  in  one  dimension 
with  Cartesian  geometry.  Executing  the  program  required 
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compiling  the  main  module  code  sphc.c  with  the  1-D  shock  tube 
setup  file  sph_init.c.  While  the  program  performed  as  expected 
for  the  shock-tube  test  case  (to  be  discussed  later),  the 
rarefaction  and  2-D  blast  wave  setup  files  did  not  successfully 
compile  with  this  version  of  SPHC.  The  code  for  these  two 
setup  files  came  from  the  SPHC  user  manual  that  accompanied 
the  main  module  code.  Failure  of  these  two  test  cases  to 
execute  with  the  available  version  of  the  SPHC  code  proved 
conclusively  that  this  version  of  the  SPHC  code  was  applicable 
only  for  performing  1-D  shock  tube  tests  and  inadequate  for 
evaluation  of  2-  and  3-D  blast  wave  problems. 

At  this  point,  obtaining  the  correct  version  of  SPHC  for 
blast  wave  simulation  was  not  pursued  because  it  required 
further  delay  of  the  effort.  Second,  and  more  importantly, 
the  setup  file  for  the  blast  wave  test  case  required  code 
modification  to  deal  with  the  configurations  of  interest  such 
as  shock  wave  stagnation/amplification  in  a  corner.  The 
expertise  was  not  readily  available  to  perform  this  task. 
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VI.  Conclusions  and  Recommendations 


A  smooth  particle  hydrodynamic  code  was  evaluated  using 
shock  tube  simulation  as  a  test  case.  A  summary  of  conclusions 
from  the  results  presented  in  this  study  and  recommendations 
for  further  investigation  are  discussed  in  the  following  two 
sections . 

6 . 1  Conclusions 

SPHC  displayed  the  capability  to  model  shocks  accurately 
at  low  shock  strengths.  This  confirms  previous  reports  in  the 
literature.  At  compression  ratios  of  interest  to  the  weapons 
effects  community  the  code  produced  an  anomalous  spike  in  the 
density  profile  at  the  contact  discontinuity.  The  literature 
has  not  reported  such  behavior  since  shock  simulations  have 
been  performed  at  low  shock  strengths.  SPHC  may  have  the 
ability  to  eliminate  this  behavior  by  use  of  a  particle 
division/recombination  method,  but  this  was  not  confirmed  in 
this  study.  The  density  profiles  generated  by  SPHC  did  not 
exhibit  the  dip  in  density  and  oscillations  at  the  tail  of  the 
rarefaction  wave  as  did  the  Lagrangian  hydrocode. 

A  drawback  of  the  code  is  that  values  of  pressure  are  not 
explicitly  calculated  as  an  output  variable.  Temperature  and 
density  values  are  calculated.  Pressure  could  therefore  be 
calculated  for  gases  obeying  the  ideal  gas  law.  Another 


drawback  is  that  the  code  is  problem  specific.  Each  version 
of  SPHC  is  dedicated  only  to  one  problem  such  as  rarefaction, 
shock  tube,  blast  wave,  etc. 

Run  times  for  SPHC  scale  between  N log  N  and  N2.  The 

Lagrangian  hydrocode  scaled  0(A/2).  Computation  times  for 
SPHC  were  four  times  greater  than  for  the  hydrocode  for  N  £  500 
particles  to  achieve  similar  resolution  and  may  actually  be 
greater  by  a  factor  of  twenty  if  the  codes  are  run  on  the  same 
computer.  However,  for  problems  involving  a  larger  number  of 
particles,  run  times  for  the  hydrocode  scale  at  a  steeper  rate 
and  will  exceed  run  times  for  SPHC. 

6 . 2  Recommendations 

Although  the  evaluation  of  the  SPH  method  was  preliminary 
in  nature,  the  technique  itself  merits  further  study  as  does 
the  MRC  SPH  code.  SPH  methods  have  only  been  in  existence  for 
a  relatively  short  ten  year  span  and  only  recently  have  they 
been  seriously  applied  to  non-astrophysics  problems.  The 
ability  to  handle  complex  flowfields  without  mesh  entanglement 
and  superior  computational  performance  in  two-  and  three- 
dimensions  make  particle  methods  a  prime  candidate  for 
multi-dimensional  blast  wave  calculations. 

Specific  recommendations  for  further  study  of  SPHC  are: 

1)  a  study  of  the  2-  and  3-D  blast  wave  and  fireball 
test  cases  performed  by  MRC, 

2)  compare  SPHC  blast  wave  test  case  to  an  existing 
blast  wave  code  using  the  same  computer  for  run  time 
comparisons , 
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3)  evaluate  the  inertial  confinement  fusion  version  of 
SPHC  for  credibility  against  recent  ICF  work,  and 

4)  investigate  the  various  physics  options  available 
for  program  execution  such  as  thermal  and  radiation 
diffusion  options. 
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